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PROGRAM DESCRIPTION GUIDE 


A. IDENTIFICATION 
Program Name 
Author 

Source Reference 

Programmer Contact 
Date of Issue 


Ordinary Differential Equation Solver - HMS 

H. Sloate, G.E. Co., Electronics Laboratory 
Syracuse, N. Y. 

An Implicit Formula for the Integration of 
Stiff Network Equations, TIS No. R70ELS-2, 
January 1970. 

V. J. Marks, GE/MSC, Houston 
May 1, 1973 


B. GENERAL DESCRIPTION 


The program is designed to provide numerical solution to systems of 
linear and nonlinear first-order ordinary differential equations. 

The program contains a new integration algorithm for the solution 
of initial value problems and is particularly efficient for solving 
differential equations having a wide range of .eigenvalues . For the 
classical methods, the integration step size is limited more by 
stability considerations than by accuracy. The new implicit fourth- 
order linear multistep method based on Gear’s formula utilized in 
this program does not become numerically unstable as the time step 
size becomes large. Since the integration formula is a multistep 
method, i-t must be started by some other means. Gill's fourth-order 
Runge Kutta method is used to calculate the three points in addition 
to the initial conditions needed to start the multistep process . 

After the starting procedure is carried out and an acceptable time 
step size found, the step size is still' kept small enough for accuracy. 

•If it is too large, it is halved. If It is too small, it can cause 
excessive running times and an appropriate test is performed to detect 
this event and the step size is doubled. 

G. • USAGE AND RESTRICTIONS 

Machine and Compiler Required - UNTVAC 1108 and FORTRAN 

Peripheral Equipment Required - Card Reader j Line Printer 

Approximate Amount of Memory - 6,027 
Requii-ed 


Qp 


v* 
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D. PARTICULAR DESCRIPTION 


Consider the system of first-order differential equations given in (1); 

^-_f (y,t), y(0) y 0 (1) 

The function^ may be a non-linear, time-varying function of^ and hence 
it may be difficult or impossible to find an analytic solution to eq. (1.) if a solu- 
tion exists. If one does exist, the above system of equations may be integrated 
point by point in time to obtain an approximation to the solution of eq. (1), ^r(t). 

There are many formulas used for the integration of ordinary differential 
equations. Some are designed to integrate differential equations of high order^ 
'while others are used for a system of first-order differential equations such 
as (1). 

The form given in eq. (1) is general because differential equations of any 
order can be rewritten as systems of first-order equations. Thus it is of 
interest to examine integration formulas which can be used to solve eq. (1). 

Classical integration techniques include the linear multi-step methods, 
the predictor-corrector methods, 2 and the Runge Kutta methods. 3 One might 
ask why new methods are being invented when we have such a repertoire of 
existing methods at our disposal. The following discussion will point out the 
shortcoming common to all classical methods. 

1 . STIFF DIFFERENTIAL EQUATIONS 

A stiff set of linear ordinary differential equations is defined as one which 
has very large and very small eigenvalues. The term stiff probably comes 
from structural engineering where stiff members gave rise to large eigenvalues 
in the differential equation formulation. 

When numerically integrating stiff differential equations, one encounters 
two conflicting requirements: a) the time step size must be large to reduce the 
number of steps taken and consequently the labor required to obtain the solu- 
tion and b) the time step must be small to prevent the integration algorithm 
from becoming numerically unstable. Stiff differential equations aggravate 
the above conflict because the maximum allowable time step for numerical 
stability is inversely proportional to the largest eigenvalue. Many time steps 
must be taken to display the solution associated with the small eigenvalues of 
the system. The combination of small step size and long running time causes 
the solution to be calculated at a great many points. 



(2) 

( 3 ) 
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Substituting eq. (3) into (2) and solving we obtain the recursion relationship: 

y . ss (1 - A h)y (4) 

J n+1 v ,J n v ' 

The solution to eq. (4) is 

y - (1 - A h) n y _ (5) 

u O 

For the sequence [Y ] to be bounded it follows that: 

| 1 - X h| <1 

0 < Ah < 2 

0<h<! (6) 

That is, if the numerical solution is to have the same property as the actual 
solution to eq. (3) (i. e. , boundedness) then there is an upper limit on h which 
is inversely proportional to the largest eigenvalue in the system of equations 
being integrated. This is generally true of all classical integration formulas. 

2- IMPLICIT INTEGRATION FORMULAS 

Much research has been done in recent years to circumvent the numerical 
instability problem. One approach is to examine the equations being integrated 
and to modify them to delete the troublesome eigenvalues. 4, 5 However, this 
may be time consuming and not completely general. Another approach is to use 
a non-linear integration formula which uses exponentials rather than polynomials 
to fit the functions being integrated. ® This type of formula has excellent stabi- 
lity properties and is exact for linear systems (no truncation error). However, 
the Jacobian of the system of differential equations must be calculated and the 
exponential matrix generated. Here we have solved the numerical instability 
problem but have replaced it with a great deal of compilation per time step, 

A tradeoff between stability and computation per time step can be made by 
using the convergence properties of the exponential power series as is done in 
the CIRCUS^ network analysis program. 


Lacking a complete solution to the problem, let us lurn our attention to linear 
multi-step integration formulas of the form: 

k-1 k-1 

y i = \ Q.y . + h ,Ly . (7) 

i = 0 i = -1 


The g's and |3’s may be determined by requiring that the formula be exact 
up to a certain order, p, if the function, y(t), being integrated is a polynomial 
in time of order p or less. If there are more a ' s and /3 's than necessary to 
make the formula exact up to and including order p,then the remaining arbi- ■ 


trary a's and /3's may be chosen to optimize other properties of the integration 
formula. ^ 


If £ 0 the integration formula is said to be an implicit formula because 
the unknown, y p appears on both sides of the equation. These equations 
have the same ftVm as the corrector formulas in the predictor- corrector 
method. If ji - 0 the integration formula, is an explicit formula, similar to 
the predictor formulas in the predictor-corrector method. 
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The implicit equations can have numerical stability for the time step size, 
h, arbitrarily large, as the following example shows. ° Consider Euler's back- 
ward formula which is an implicit formula of order one and step number k - 1. 


n+1 


= y + hy 


n-t 1 


( 8 ) 


If we use eq. (8) to solve eq. (3) we obtain the difference equation: 


y n,l = ^ + 


Ah) 


-1 


( 9 ) 


The solution to eq. (9) is 


y = ( 1 + A n) 
n 


( 10 ) 


If Re j A j’S 0, then ,\ y ] is a bounded sequence for all h>0. Thus there is no 
upper limit on time sPep size due to numerical instability problems. Here is 
one formula, at least, which has the desired numerical stability properties. 

The question arises as to whether there are more formulas with the above 
stability properties. Before answering this question let us present the concept 
of A-stability. ® A k step method is called A-stable if all solutions of eq. (7) 
tend to zero, as n->cc, when the method is applied with fixed positive h to eq. (3) 
where A is a complex constant with positive real part. A-stability requires 
that if the solution to the differential equation is stable then the corresponding 
solution to the difference equation be stable also. The example given in eq. (8) 
is an A-stable method. Unfortunately there are not too many more of them. 
Dahlquist® proved two theorems important to our discussion: 

Th, 1, An explicit linear multistep method cannot be A-stable. 

Th. 2. The order, p, of an A-stable linear multistep method cannot 
exceed 2. 


Theorem 2 indicates that if we require A- stability we will be able to use only 
formulas of low order and hence limited accuracy. If the restriction of A-stability 
is removed, we can find high-order formulas with acceptable numerical stability 
properties. 

In Appendix I it is shown that if we want an integration formula which remains 
stable as h-»-cc, we must use an implicit formula and its maximum order cannot 
exceed its step number k. Such an integration formula will not be A-stable, how- 
ever, if p > 2. 

3. INTEGRATION FORMULA SELECTION 


As h — »-oo it would be good to have the roots of the integration formula tend 
to zero so that the unimportant large eigenvalue modes of a system of differential 
equations would rapidly die out as the time step size is increased. This is the 
approach taken by Gear 19 t n obtaining implicit integration formulas of order 2 
through 6. In eq. (7) there are 2k+l arbitrary constants. If a kth order fit is 

k j 

required, k+1 of the constants are fixed. If the k roots of v 0. C = 0 are 

i - Q 1 

all to be zero, the remaining k constants of eq. (7) are fixed and the integration 
formula is unique. The formulas of order 2 through 6 are given in Appendix II. 
Generally speaking it is desirable to use as high an order formula as possible 
to achieve greater accuracy. The higher order formulas have the disadvantage 
in dial more past values must tie stored and more manipulation is required when 
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the stop size is halved or doubled. A fourth-order formula was chosen for 
integrating the systems of differential equations as a compromise between 
accuracy and complexity. It is different from Gear's fourth-order formula 
in that its coefficients were selected to achieve a compromise between trunca- 
tion error and roots close. to zero as h-*-oo. This w^s done by minimizing a 
performance index using the Fletcher Powell method. ^ The performance 
index was a weighted sum of the truncation error given in eq. (12) and the square 
of the sum of the squares of the roots of eq, (7) as h— p-oo. 

4 . TRUNCATION ERROR 

The truncation error of the integration formula in eq. (7) is due to the fact 
that the formula has only a finite number of terms. This error would be present 
even if there were no roundoff errors in the computer due to using numbers with 
a finite number of bits. An expression for the truncation error can be obtained 
by integration by parts. 11 (See Appendix III.) It is given by 

(n+l)h 

(k+l)/ N \ f 

T n = JL Ti ~ J G(s)dS (11) 

(n-k+l)h 


where G(s) is called the influence function and is of the same sign over the 
interval from (n+l-k)h to (n+l)h. The k+1 th derivative is evaluated at N which 
is somewhere in the interval from (n-k+l)h to (n+l)h. If G(s) changes sign over 
the interval, the First Mean Value Theorem *2 which was used in obtaining eq. (11) 
does not apply and an estimate of the error can be obtained from: 


_ .y^W 

n IT! 


/ 


(n+l)h 


>L 


G(s)| ds 


( 12 ) 


(n-k+ l)h 


Equations (11) and (12) were used to obtain the truncation error of the inte- 
gration formulas studied. 

If G(s) is of the same sign over the interval and eq. (11) applies, an alternate 
shorter method can be u9ed to compute T^. Assume the error is given by 


T n " lEIfT hk+1 V <k+1) (N) < lla > 

Assume Y(t) = l and substitute into eq. (7). Since eq. (Ha) is the error in 
eq. (7) E^ can be evaluated. This is the standard way to find error terms in 
quadrature formulas. Note, however, that it is valid only if G(s) is of the same 
sign over the interval covered by the integration formula. (See Appendix III 
for an example of the procedure for calculating E]£. ) 

5. STABILITY BOUNDARIES 

Let y = qy . The roots of the difference equation given in eq. (7) are the 
zeros of 

k-1 

v (a ^ I qh;ij) ^ 1 ’ - 0, a 

i - -1 


-1 


(13) 
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Solving for qh gives 


qh = 


k-1 
- V 


k-1-! 


a .K 


i .. -1 

FT~ 


S' 

Li 




k-l-i 


i = -1 


(14) 


If there is some value of qh which causes £ to have unit magnitude, then, 
we can let t = e and let 0 vary from 0 to 7r . This will map the corresponding 
stability boundary in the qh plane. If the stability' boundary lies wholly in the 
plane HE[qh] > 0, then the integration formula being considered will be A-stable. 
If not, then there are some values of qh in the plane Re[qh]<0 for which the 
roots are not less than one in magnitude and by definition the method is not 
A-stable. 


6 . COMPARISON OF TWO FOURTH-ORDER FORMULAS 

Gear's fourth-order formula has a truncation error given by: 


<«> 

The truncation error of the optimized fourth-order formula given in 
Appendix H is: 

| T n| =1 3T^ h5 | yC5)(7? )l (16) 

These expressions can be obtained by integrating the influence functions shown 
in Figure 1. See Appendix HI. 

The stability boundaries in the qh plane of are shown in Figures 2 and 3. 
Both methods are A {a ) stable in the sense of Widlund. 13 Gear's method has 
an o -12 degrees. The optimized method has an a =63 degrees. It can be 
seen that there is a trade-off between truncation error and the size of the 
stability sector in the qh plane. The optimized formula is more accurate than 
Gear's method but its stability sector is smaller. 

7 . SOLUTION OF THE IMPLICIT EQUATION 

The integration formula given in eq. (7) is implicit if (3 ^ ^ 0. Since the 

unknown, y appears on both sides of the equal sign, the equation will have 

to be solved 1 iteratively. One approach which is used in the predictor-corrector 

methods is the method of successive substitution (Picard's Method). The k th 

estimate of y is used on the right side of eq. (7) to evaluate y n+ ], , and this is 

used to evaluate the k+1 th estimate of v , 

■ n+1. 


-1 


h y 


00 

m 1 


k-1 

V 

u 

i- 0 


( Q ; 


y . + 

n-i 


h/3. y 
l n- 


(k-t-1 ) , 

Vi ' '* 


(17) 
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If .y^°' is the initial guess for y n+ j, the solution to 'eq. (18) is: 

This is a divergent sequence if jqhj >1. This is the same restriction we 
were trying to avoid by using implicit methods in the first place. It can be 
shown that this restriction holds for systems of equations and nonlinear equa- 
tions as Well. ® 

The Newton Raphson method can be used to solve implicit equations such as 
eq. (7). It does not limit the size of qh. If eq. (7) is applied to systems of equa- 
tions such as: 

S ^n+1^ ~ im+1 ~ ^-1 h i2n+l 


. I 


[a. y . + h fi. y .) 
v l^- n-i l^-n-r 


it can be seen that g (y .) will be zero at the solution. If g is expanded in 
a Taylor series and onlyrne first term is kept, we have 

Q (y + Ay .) = Q (y^ k \) + J • a y , 

2 iin+l ~n+r w n+1' ** • -s-n+1 

Solving eq. (21) for ^ n+1 gives 

Ay . = - J 1 Q (y (22) 

~n+l ^ Xiin+r 

where J is the Jacobian of g, The Jacobian can be fo ind from eq. (20) and is 
J^-^jhA (23) 


where A is given by 


r.l y, = f 
3y i ’ y = £ 

l u y 2 . 


y = f 

n n 


' If the differential equations being integrated are linear, the matrix A does 
not change during the problem and needs to be calculated only once. For very 
nonlinear equations A may change rapidly and a new Jacobian may have to be 
calculated at each time step. In practice, however, a new Jacobian is not cal- 
culated this often. If the number of iterations required for the Newton method 
to converge exceeds a certain number, the Jacobian is recalculated using differ- 
ence techniques, it is inverted, and a new Ay n+ j is found. If the number of 
iterations per step does not exceed the limit, the inverse Jacobian from the 
previous step is used. In this manner considerable time is saved which would 
be used in evaluating the derivatives jhtl- 
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If after calculating a new Jacobian the Newton method still fails to converge, 
the step size'h is halved and eq. (7) is applied once more. If after b has been 
halved a certain number of times there is still no convergence, the simulation is 
terminated and an error message printed out. 

If the system of differential equations being integrated is very stiff and the 
transient has died out, qh»l for all the eigenvalues of the system. 

[I- A] - j-yTr A' 1 . (24) ‘ 

Equation (24) shows that inverting the Jacobian - is almost the same as inverting 
A which has extreme ranges in eigenvalues. Due to roundoff in the computer, 

A seems nearly singular and give's very inaccurate results when it is inverted. 
An iterative scheme to correct the- elements of the inverse matrix has been used 
to overcome the problem. ^ It greatly reduces the number of iterations used 
for convergence of Newton's method in some cases. 

8. SIMULATION PROGRAM ORGANIZATION AND LOGIC 

a.. Starting Procedure 


Since the integration formula in eq. (7) is a multistep method it must 
be started by some other means. Gill's-^ fourth-order Runge Kutta method was 
chosen to calculate the three points in addition to the initial conditions needed 
to start the multistep method. Two more points are calculated using the multi- 
step method and these six points are used to approximate the fifth derivative of 
y(t) by difference techniques. From this the truncation error is calculated and 
the accuracy test is made to determine if the time step size is too large. If it 
is, the step size is reduced by a factor of ten and the starting procedure is re- 
initialized. The error test is: 

| T n| <(1+ IVl| ) * ELIM <25) 

If this inequality is satisfied, the time step size is acceptable. ELIM 
is data input by the user. If y(t) is an extremely large number, roundoff error 
in the computer will keep T n from being less than ELIM. The ly J term is 
included to prevent the occurrence of this situation. I i 

b. Time Step Size Control 


After the starting procedure has been successfully carried out and an 
acceptable time step size, h, found, eq. (25) is still used to keep h small enough 
for accuracy. If it is too large, it is halved. If it is too small, it can cause 
excessive running times. Hence the test in eq. (26) is made. 

[T n | < (1 + | y n+J | ) * ELIM/40. (26) 

If this inequality is satisfied, the step size is doubled. The factor of 
40. was chosen because the error is proportional to h . Doubling h multiplies 
the truncation error by 32 if y^(/y) remains constant. The factor of 40 was chosen 
so that eq. (25) would still be satisfied after h was doubled. 
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Since six points are necessary to calculate the fifth derivative, it is 
necessary to have five accurate past points before y n - h i is calculated. Before 
doubling h,eight past values of y(t) are required so that there will be five past 
values of y(t) when y n+ j is calculated on the next step. (See Figure 4. ) The 
step-size doubling is inhibited until eight past points have been calculated. 

If Yn+l is calculated and the truncation error is too large, the step size 
is halved. (See Figure 4. ) yn-1 and y n „3 are found by interpolation, y n _^ and 
y „ are calculated from these values using the state equations. The old y n -l, 
yn-1 are placed in y h _ 2 , y n -2. The old y n _ 2 , y n -2 * re placed in y n _ 4 ,, y n _ 4 > 

A new yn+l is calculated from the integration formula. If the truncation error 
is still too large, the step size is halved again and the process repeated. 


n-4 n-3 


n-2 n-1 


n 


n+1 


X X 

X X X X 

n-7 n-6 n-5 n-4 


XXX 

^ double step size 

X X X . X X 


n-3 n-2 n-1 n n+1 


X next step 


DOUBLING THE STEP SIZE 


n-4 n-3 n-2 n-1 

X X X ' X 


n n+1 

X X inaccurate 

point 


halve step size 


XXX XX 

n-4 n-3 n-2 n-1 n 


HALVING THE STEP SIZE 
Figure 4. Logic for Changing Step Size 


c- The Newton Iteration 


The Newton Raphson method is used to drive Q to zero. An initial 
estimate of the solution, yn+i, is calculated by extrapolation using the five 
previous points which have been stored. Using this estimate §1°' is calculated 
and Ayn-i-l is calculated. 


(k+i) 00 

n-»l n fl 


v TI 


Ay 

^n+l 


(27) 
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The new value of^Vl is computed using eq„(27) with TI = 1. is 

evaluated and if ^ Q^°\ a new Ay n+ |is calculated and the process 

is repeated. If not, THs decreased oy a factor oi 10 and eq. (27) is applied again 
with the same Aj^i+1- This process is repeated up to NDEX times. If there is 
still no success^ a new Jacobian is calculated and a new Ay n+ -j is used in eq. (27). 

If there is no success even then, h is halved and the whole^ procedure is tried 
again, h can be halved up to NHLIM times. If no success is obtained, an error 
message is printed and the program terminated. 

New A^ n+ i values are checked in the inequality 

|Ay| < (SC + |y|) * XNLIM - (28) 

If eq. (28) is satisfied for all the state variables, the iteration has been 
completed successfully and the last estimate of v n+ i is taken to be the solution 
of the implicit integration formula. SC is usually set to 1. and XNLIM is usually 
set to 0. 1 * ELIM. 

Each time a Ay n+ i step has been successfully carried out, the estimate 
of the inverse Jacobian is updated using Broyden’s scheme. 1® This updating 
keeps the estimate of the inverse Jacobian current and lessens the need for 
calculating new Jacobians. 

F. description of input 


The differential equation solver is a stand-alone program. Flow charts of 
the program are shown in Appendix IV. It needs two sets of input from the user . 
The first is input data which is read in on three cards. 

The first card has a format (2F8. 1, 15) and reads the variables TF, PT, 

NV. TF is the final time for the problem. When the solution time reaches this 
value, execution will terminate. PT is the number of points in time the user 
wants printed out. For example, if TF = 10. and PT = 5., the state variables 
would be printed out at t = 2. , 4. , 6. , 8. , and 10, seconds. NV is the number 
of state variables. 

The second data card has a format (10F8. 2). The initial values of the state 
variables (up to 10) are input here. If there are less than 10 state variables, 
leave the unused columns blank. 


The third card has a format (4E7. 1) and reads the variables YLIM, ELIM, 
XNLIM, SC. YLIM is usually set at about 10~ 5 and governs how the state varia- 
bles are perturbed to calculate the elements in the Jacobian matrix. The program 
uses YLIM as follows: 



< YLIM, perturb y^ by 10 .. 

-3 

< YLIM, perturb y^ by 10 *y u 


ELIM controls the truncation error of the integration formula. See eqs. (25) 
and (26). It is usually set at about 10~ 3 . If a more accurate solution is required, 
it can be set smaller. However, a smaller ELIM requires more time steps to 
complete a simulation. XNLIM controls the convergence requirements for the 
Newton iteration. See eq.(28). It is usually set to 0. 1 * ELIM. 
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The other input required of the user is a subroutine in FORTRAN which 
describes the system of equations being integrated. An example is shown 
below. The maximum number of state variables is 10. The state variables 
are contained in the array X. The first derivatives of X are stored in the array D. 
The variable NUM is used to count how many times the state equations have been 
evaluated during a simulation. 


1 

2 

3 

4 

5 

6 

7 

8 
9 

1 0 


CSTATE 

SUPROUTJNfc STiT£(X,[n 
COMMON Y0»|1Y(V.DB.T,H,NV,NUK 

double precision yocioi.dydcio) , l» b c i n j 

DOUBLE PRECISION X<i),D<l> 
D(2)=X(1)-X(?) 
utl)**in()fl.*(x(i) ♦ p ( ?) > 

Ny M = NyM + j 
RETURN 

end 


F. DESCRIPTION OF OUTPUT 

Since the integration formula is a variable step method, the variables 
may not be calculated at the value^of time at which printout is desired. Hence, 
it is necessary to use Lagrangian interpolation to find the values of the state 
variables at the desired time. See Appendix B for Sample, Output. 

G. INTERNAL CHECKS AND EXITS 
None. 

H. INDEPENDENT SUBROUTINES 
None 

I. SYSTEM SUBROUTINES 
No special subroutine. 

3 . COMPLETION OR FINAL CHECKOUT DATE 
April 25, 1973 

Source listing and sample cases are given in Appendix B. 
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MATHEMATICAL PROOFS 

DEFINITION - A linear multistep method is said to be consistent if it is of 
at least order one. 

k i 

If a linear multistep method is consistent, £ a. £ has a simple root at 

i-0. 1 

£ - 1. Let y •• K, a constant. y = 0. Sustituting this into 
k k 

2 V y n + l“ h £ ^i y n+i * 0 • (1) 

i=0 i=0 

results m 

£ «, = 0 <*> 
j.O 1 

If there are two roots at £ = 1, the solution of the difference equation is 

L- K i + K 2 n + ¥s” + - + \ c (3) 

where the K's are determined by initial conditions. The sequence [y ] is not 
convergent as n -»co because of the term containing n. Hence the metnod will 

k i 

not be convergent if ^ q ( has more than one zero at £ = 1. (See Henrici, P. 

i=0 1 

Discrete Variable Methods in Ordinary Differential Equations, Wiley, New York, 
1962, pp. 217-218.) 

k 

We will now show that consistency requires that ^ P ■ 0. Let y = 1, 


y(o) = 0. The solution to the differential equation is y(t) = t . Equation (!) becomes 
k k 

Yj q i • i " Z Pt = 0 ( 4 > 

z r i 1 ,• * 


. k 

Since V a £ 1 

d-C i-0 1 


k 

eq. (4) that \, /i. / 0. 
i‘ 0 1 


j? 0 because £ = 1 is a simple root, it follows from 



Lemma 1. If the following three hypotheses are true: 

1) £ — - Z », z 1 

£ b, Zl ‘-° 

i=0 1 

2) b = 1, b >0 i = 1,2, .... 

2 

3) b , b , - b >0 
n+1 n-1 n 

then a = 1 and and a.<0 i = 1, 2. . . 

o i 5 

Proof: 


CC 

V 


CC 

S’ 


If ' . „ has a non-zero radius of convergence, then 

i=0 i Z • ' i = 0 


a i Z 


(See Knopp, K. Infinite Sequences and Series, Dover, 1956, p. 116. 


From 


CC CC 

^ b. Z 1 . V a. Z 1 = 1 

i=0 1 i=0 1 


we obtain the following: 


b a = a =1 
oo o 

b l +a l = 0 

b 2 +b l a l + a 2 = 0 


b + b .a, +. .. +b, a . + a =0 
n n-1 1 1 n-1 n 


Solving for b^: 


n 


b =~^ b . a. 
n . , n-j ] 
j = l J J 


In a similar manner 


^n+l .t-L b , . a - a . 
j = l n+l-j j n+1 


does also. 

) 


(5) 


( 6 )’ 


Multiplying eq. (6) by b n and eq. (5) by b n+ j and subtracting gives 



Since by hypothesis - b • >0 


b b. b- b . b 

b, b, ^b„ ^ "d . . 'b . 

12 3 n+1-] n+1 

b . b . - b , b > 0 and this shows the term in paren- 
■ n-j n+1 n+1-] n 

thesis in eq.‘(7) is always positive. Equation (7) will now be used to show by 
induction that a j < 0, i - 1, 2, . . . 

a ' -1 
o 


from eq. (7) 


a 1 a - b 1 <0 


b l a 2‘ !b o b 2 ' b l )a l <0 


Assume 


a 2<° 


a r a 2 ] 


, a^< 0. From eq. (7) 


1 n ' 

a 1 =-j- — X (b • b . - b , . b ) a. 
n+1 b n ^ n-j. n+1 n+l-j n' ] 


J 

Since all the terms in parentheses are positive, b n is positive, and the a. 1 
are negative', a n +l must be negative. This completes the induction. J 

Theorem : No explicit linear multistep method remains stable as h becomes 
arbitrarily large, p = k is the highest order possible for an implicit linear 
multistep method which remains stable as h becomes arbitrarily large. 


Proof: 


let y = qh 


( 8 ) 


The linear mulhstep method used to solve eq. (8) can be written as: 

0(hPtl, 


1=0 


Substituting the solution y = of eq. (8) into eq. (9): 


k 

V 

L-t 


i=o L 1 


a. (e qh ) n+1 - qh a (e ph ) n+i 


= o<h p+1 ) 


(9) 


(10) 


qh 


Let us use the transformation ^ =■ e H which maps the LHP of qh into the 
unit circle in the C plane. Dividing eq. (10) by qh and substituting C = ed* 1 and 
qh - InC: k 

k 


y ... 


• 1-0 


1 'l-t 


1 ! 1+Z 1 
i„ ( T - f -) 


- v /) t; 

i-0 


O (| int'l p ) 


( 11 ) 
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The C 's may be viewed as the roots of eq. (9). For stability they must 

r_ i 

remain in the unit circle. The transformation Z —maps the interior of the 

C+ 1 

unit circle into the LHP of the Z plane. The inverse relationship is 

< (12) 

Substituting eq. (12) into eq. (11) 

.^. a i tz* Vo ,1+Z ^ , n ( 1+Z« P , /. o \ 

i — =? — 4 i-z ' 0 ( f ln ( 1-Z ^ ^ (13) 

t f l+Z, i=0 
In (j^-) 

1-Z K 

As h-*0, C *■►1, and Z-^0. Therefore multiplying eq. (13) by (-^ — ) and 
simplifying: 

■ -- ( Z) . S (Z) = 0 (Z P ) (14) 

W-ra-) 

where 

r(Z) - (I#)" l a, £* )‘ . I a. z‘ (15) 

S(Z) = l p &*) = £ b z 1 (16) 

1 i=0 ll-Z i=0 1 

If eq. (T4) is to be of order p then the coefficients on the left side of eq. (14) 

p v i 

must cancel up to the Z term. r(Z) has a zero at Z = Obecause a. K has a 

k i=Q i 1 

root at C = 1. Since eq. (9) is stable for h = 0 all the roots of £ a. £ must be 

i=0 1 

on or within the unit circle and hence a. > 0 in eq. (15). If we desire that all 
the roots of eq. (9) be on or within the Unit circle as h-*-oo then b. > 0 also. Let 

— Tf g- = c 0 + c 2 z 2 + c 4 z 4 + ... (17) 


Due’ to Lemma 1: 

C 0 >0 > C 2» <0 > V = l > 2 > 
Matching coefficients in eq. (14). gives 

b 0 = C 0 a l 

b l C 0 a 2 


( 18 ) 
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b 2 v = C 0 a 2 v+1 + C 2 a 2 u-1 + • ‘ * + C 2 i/ a l 
b 2(/+l = C 0 a 2 u+2 + C 2 a 2 i' + ' * ' + C 2^ a 2 

Since S(Z) is a polynomial of kth order b n - 0, n > k. The same applies 

fur a 0, n > k. 
n 

Assume k is odd. 


2 v+l - b k = c Q a k+1 + c 2 a fc _ 1 - + . . . + c k a 2 

Since ak t-1 - 0, the coefficient matching in eq. (14) requires that b,< 0. But by 
Descarte’s Rule of Signs this means that S(Z) has roots in the right half of the 
Z plane and hence f C|> 1 as h-^ oo. The conclusion is that we can only require 
coefficient matching in eq. (14) up to the k - 1 th term. Hence p = k. The same 
reasoning holds if k is even. 

k 

If eq. (9) is explicit, then 0, = 0. From eq. (16) S(l) = 0 = £ b. . Since 

K i=0 1 

bj > 0, i - 0, . . . , k if all the roots are to be within the unit circle as h-r <», 

* jk k 

all the bj’s must be zero, including b Q . But bQ=(?£) £ 0^ from eq. (16). 

i=0 

From the consistency condition 

ifo 

1=0 1 

Hence if 0k = 0 and if we require stability as h ->qo, we find the consistency 
condition is violated. The conclusion is that no explicit linear multistep method 
can be stable as h -*■«>. 
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GEAR'S FORMULA 


k -- 2 


. k = 3 


k = 4 


k = 5 


y n + l =!*„- 5 Vl +h <!’W 

y n+l 11 y n ~11 y n-l + TT y n-2 + h ( 11 y n+l^ 

„ _ 48 36 16 

y n+l 25 y n ~ 25” y n~l + 25 y n-2 

‘ A V3 +h ( S y n + l } 

_ 300 300 200 

y n+l = 13 ? y n " 137 Vl + 137 y n-2 

“13T V-3 137 y n-4 + h ( 137 y n+l^ 


k = 6 


n+1 


= 360 
147 y n 


450 400 

147 y n-l + 147 y n-2 


- 225 72 10 

147 y n-3 + 147 y n-4 ~147 y n-5 


+ h ( 


60 

147 



) 


Coefficients of Optimized Formula (k = 4): 


o Q = 1. 584 

*-1 

- 0.4539 

o 1 = -1.017 


= 0. 235 

a 2 = 0. 529 


- 0.0568 

= -0. 0968 

@2 

= 0. 00567 


P 3 ~ 0 , 000201 
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TRUNCATION ERROR 


An error term for a Taylor series expansion of y(t) about a point in time 
at nh can be obtained from the integral 

(n-i)h 

p- f ([n-i] h-s) r y^ r+1 \(s) ds 


nh 


If this is integrated by parts, we obtain 
(n-i)h 

FT / [n-i]h-s) r y <rtl) (s) ds =4. ([n .i| h . s) r y (r) (s) 


nh 


r : 


(n-i)h 

nh 


+ 1 (n-i)h r _i 

7^-jyf / (|n-i]h-s) 1 y { ' (s) ds 

v ' nh 


(n-i)h 


y^ (nh) + icr^Ty*" / ([ n-i]h-s) r_1 y^(s) ds 


nh 


If integration by parts is carried out r times and the terms rearranged: 

2 2 

y([ n-i]h) = y(nh) - ih y^ (nh) +-^r — y^ (nh) - . . . 


OJW y w (nh)+ _( /n ;‘ h . s) y^> (s)ds (1 , 


(n-i)h 


A Taylor series expansion with error term can be obtained in the same 
manner for the first derivative of y(t) about a point nh. 

y (1) (| n-i | h) = y (1) (nh) - ih y (2) (nh) + 


, ..r-l.r-1. r-1 , . 

. 1-0 i h Jr) 


i (n-i)h 

( r _^y; y lr) (nh) + -( pn)7 /„h (J n-i]h-s) r "y tr+1) (s)ds (2) 


Equations (1) and (2) can be substituted into the integration formula in eq. (3) 
and an expression for can be obtained. 


k-1 
_ \' 


, - a ■ y . + h 

n+1 u i J n-i . — - 

l-O i=-l 


k-1 

T p.y . + T 
- 1 f l J n-i n 


( 3 ) 


Since eq. (3) is exact up to and including order p we obtain for T^: 
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T 


n 



( f n+ 1 ] h-s) p - ph/3_ j (| nf f]h-s") r * 


k -1 



(( n-i ] h-s) P + phfi j ([ n-i jh-s) r_1 ]| 


y^ p+ 1 ^(s) ds 


where 


(4) 


([n-i] h-s) -i 


[ n- i j h-s 
0 


[ n-i ]h£ s< nh, i ^-1 

nh < s, i = -1 
otherwise 


Equation (4) can be written 


- /*(n+l)h . 

T n = — j- / G(s) y tP+ 1 ) (s) ds (5) 

-'(n+l-k)h 

If G(s) is of the same sign over the interval £(n + 1 - k) h, (n+l)hj, then 
J G(s) ds is a monotonic function. If yf^-^s) is continuous over the interval, 
then the First Mean Value Theorem applies. * 2 The estimate for the truncation 
error can be written: 


(n+l)h 

T n = y P+1 (V)[ G(s) ds ( 6 ) 

*fn+l-k)h 

where 77 e [ (n+'l-k)h, (n+l)h ] . 

Example 


Let k = p - 4 

y n + l = °o y n + °l y n-l + ff 2 y n-2 + G 3 y n-3 

+ h ! ' 3 -l4,l * >Vn * |! lVl + 'Vn-2 + 'Vn-sj + T n 

The influence function is obtained from eq. (4). 

G (s) = Gj (s) = a 3 (| n-3 j h-s} 4 + 4h/? 3 ({ n-3 ] h-s) 


(7) 


for j n-3|h< s < | n-2 ]h 
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G (s) G 2 (s) - G t («) + o- 2 (| n-2 |h-s) 4 t 4h/i 2 (.( n-2 |h-s) 3 
for | n-2 ] h < s < [ n-1 | h 

G (s) •= Gg (s) = Gg (s) 4 c ^ ( [ n-1 ]h-s) 4 + 4h/3j ( [ n-1 ] h-s) 3 
for ( n-1 ]h <s <nh 

G (s) = ([n4l ]h-s) 4 - 4h/3 ^ ([n+l]h-s)^ 


for nh< s <[n+l]h 


For Gear's Method 
48 

a_ = -STT , o, = - 


.0 " 25 ’ 1 


36 

25 


a n = 


16 

25 


= 


25’ ^-1 " 25 


12 


'*0 “ ^ " *2 " *3 = ° 

The influence function is shown in Figure 1. Since the influence function 
is of the same sign over the interval ([n-3]h } [n+l],h), the error is 

T „ “ IT ' 1,5 

5 

To find E 4 assume y(t) = t . Substituting this into the integration formula: 

5 .5 „. 5 , 


5 „ ,5 

-c 

5 


h v = -a x h w -32a 2 h -243<* 3 h + 5hO _ 1 + ^ + 16/3 2 + 81^) 


+ E 4 h 


Solving for E^: 

E 4 = 1 4 a x 4 32a 2 4 243« 3 -5^^^+ 16^ 4 810 g ) 
Using the values of the a's and fi's for Gear's Method gives 


E 4 = -11.52 


2. 304 .5 (5). , 

= - 51 — h y (v) 


This shows that 


/ 
\J( r 


(n + l)h 


G (s) ds 


= 2. 304h' 


4(n-3)h 

which can bo verified by integrating the influence function directly. 
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; 1 •' ■ C.M A I N - • 

— 2* — s — ' — C 0 -H O-y-0- r D B t ^ V f Mli-M r Y | --Y-1 rV 2 -y V3 — f€-j- ) DY >- P Y t t : D Y2 t 0*3 j - 

3* ' ■• 1 DYV t TPR,DT,NDEX f> IACiO,G,N5 ' ' .' • " 

M-*- .......... CO-M-MO Hy B N D -5-/-Y-!-!. M»-.K -N LT— •£ L - 1 - M | Sc.tMEQ3 i N H — — — - — — — — - — — ■ — .. — ~ 

•s'* DOUBLE PRECISION YOIIO) ,OYOUO) »DBt!0) *YUO) »Y1 f 10> t V2U0) »V3( 10) 

6* J- - yil-tt )- i-O-Y- (10J * 0 Y- 1-C 1-0 )- # O-Y-2 LMU-* oY-3-(J 0J-*4)^X4-C-)-,£)-U-0-L4.&-(-J-O-»-i-OJ-,.- 

• 7* 2 ‘ Y5'( JO) >Y6< 10) ,Y7( 10) ,DYSC 10) iDYfr! JO) »DY7( 10) 

- — a-*- --- ~~0 O.u B L-E--P-R E-CT-5-I-0 3 MA X ,Se-T : 

■9 * ‘ DATA LJM/9/»HM{N/| tOE-15^ 

— ~1tQ-*~ C — — : — - 1 — — • 

n* C I Ml T | A u I ZE - NF COUNTS STARTING FAILURES 

-12* - r C 1 naL-C-OUN-T-S — 5-UGC-E-SSFU-L — I-N-TE-G-R-A-T-J-ON-S T — 

.13* C N6 COUNTS S-TEPS SINCE- LAST HALVING OR D0U8L I NG 

. C N H- F-u-ASS— C-oN-V ER-GE-NG-E — EN— F-P^ - — 

• ■i-5* c.' lih-is’the starting failure limit. 

— T 6*-- G- ; — : v cM kS -.- C -O- U -N-T-S-H a - T R -I X ■ I NVERSIONS ■ JAC O B I A N -£ V - ft L -- U A - TTO fH 

1 7 *‘ C ■ NUN COUNTS EVALUATIONS OF STATE EQUATIONS 

■ IS* ' • 1 PR-lVfT'~3 J *■ ----------- — : : 

19* 3 FORM AT ( 2 | Hi I N.I T J A L C OND I T I ONS ) 

■ -2 0* : ; -A-S-S-TG#— 9 8-T-O-NSW ' 1 ■ 

2j* NH » 0 

—2-2-* — ; NF — m — 0 ; 

23* JAC * O' 

-2M* t NUN-« -Q : ' : 

25* . READ S,tf,PT»NV 


-2-6-* "-S--F l 0R-M-AT-(2F8 »-^T-5-> 

27* NEQ3 «. MAX0(NV*3.7) 

-2-0* N-D-E.-X--« — 6 — r— 

29* DT 5TF/P T 


3 J * . H * TF « o *001 

- 32-* . R E-AlD — 6-^J-Y. 0 -L I-> -»4-*-l. t -N.V-} 1 . 1 : 

3.3* 6 FORMAT! 10F8.2) 

- • -3 -0* P-R-l-N- T - — 0 t T - Y - O - f-T -) — ) - rW -V-) — 

.35* H FORMAT! I X , 6 E l 5 - -7 » 

3 A * R E-A (5 — L0~ Lt'Y~LT M-y E-L T M i X ^ LT M-f 5 C ■— — — — 

37* J 01‘ FORMAT l ME 7, 1 ) 

--3 -g-*- — — — - — 1 P-R'J-NT- — f-0 0- * * — : — - — 

39* J DO FORMATJ9H0 LIMITS) 

-"^0* r-F-ft-i- NT- -M r — -*4rT-H : , CLi H ; 

Ml*- • E L I M » EL ! M * SC 

- M 2 *- E-DUB «--EtlMy- , f(y. 

M3*. CALL-STATE ( Y q » D Y 3 ) 

-M*)-* — - - - -.do -26-1 -N-V- : 

MS* 26-1 T3 ( I) * Y o < I ) 

— M-6* — € — : — —— — 

M 7 * C U 5E GILL’S METHOD TO -START 

- -9b*- - — 10 — T — — Q-r- - — • — — - — — * — ■ — - — • — - — ■ ■ ~ . — — - — - . -■ — — — — — ■ - — — — — - 

' 99* ' * ‘3 . 

.. /.So-* N-6--S— 3 ; : ; 

5 : if CALL GILL ’ . 

-..-5.2-4- ^ 00—2^2 — N-V 1 

- S3* - - pY2( I ) * dYq( I > 

- &N*- 262 - -Y-2-I -Y-O/C I -) : — — ‘ 

'55* . CALL GILL . . 

-56*- -00-263- !-»-l-»-NV- ^ 

- S’ 7 * . D Y 1 ( I ) * o Y 0 < I ) ’ 

8-> - -26-3- — V-l-r-f— ) — »~¥!34rh4 — — — 

• 5-9 * "call gill • - ' 1 
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6i* 

— 6 - 2 -*- 

63* 

- 6H-* 

65 * 

-. 66 * -- 
67 * 

— 6 . 8 -* 

69 * 

? 0 * 

71* 

72 * 

73* 

- 7-J) * 

78* 

-76* - 
77* . 
78* - - 
7 9 * 

~- 8 Cl * 

8 1 * 

. 82 * • 

93* 

8ft* -- 
88* 

— 8 6 * ■ 
87 * 

. SB* - - 
89* 
9-0-*-' - 

? l * 

. - 9 - 2 -* 

93* 

-9 4)*. . . 

95* 

. 9 . 6 * . . 
97* 

--9 ft-* 

99* 
1 . 00 * - 
101 * 
10-2* - 
103* 

10 - 4 * 

1 0 8 * 
106 * * 

1 07* 

1 OS'* ■ . 
109* 

t 1-0-*-— 
111* 

I 1 2* 

II 3*. 
Ill* 

1 1-5* 

1 1 6-* — 
117* 

11- 8-*- 


C 

■ C. 


INITIALLY, the 8 POINT SCHEME IS USED 

00 — i _ 10 — 1“ *)~r NV — ■ ■ — ~ — 

llo YU) - ft » * < Y 0 ( I ) * Y2 ( I ) J -6,*YI ( I )-Y3tI) 

- T--«- T-+ H 

C-A-L-C-UL A-T-E -I-N+X-l-A-L- -JAC-OB-W N- 

CALL STATE ( Y , 0 Y ) 

CA LL— Q-X (-OB-.-Q-M-A-X-) 

00 13 J * 1 , NV 

X * Y { 1 ) 

DIE * DA8SU) - YLIM 
IF (DIF) l6,l6rt-5 • 

16 YU) ■ X ♦ 1 • 00-ft 

S E~T — m . _( *-004 

GO TO lft 

15 Y.U ) « .X*-l .00-1 00 

SET * 1000*/X 

j.8 CAL-U-STA-TE ( Y-, 0 Y ) - 

CALL QX(Q,QMAX) 

~ 0 . 0 — )~ft . - j .«. l -^. y - - 

1 ft G ( J , I ) a { Q ( J ) - OB ( J ) ) * 5 ET 

-4.3 .UU.t.X - 

call matinv <nv,G) 

• J AC- -JAC- +-1- - 

GO JO 107 


-- C - 

c 


c 

c 


-2 2 2 


H- j 

N 6 


TO HERE 
1 . H-+H- - - 

a ft 


WHEN H is TO Be doubled 


SHIFT Y * S - 
D - Q — 1 - 60 — I - a - l - t - N - V - 


Y l is ALREADY IN POSITION 


dyoi I ) 

- D-Y 2 M ) 

D Y 3 ( I ) 1 
. P.Yil.LlJ— J 
Y0( I) * 


DY (I ) - 
- D Y-3-U-) 
DY5 t I ) 
-D.Y.7.I U. 
Y ( I ) 

-Y-3-UJ- — 


Y 3 U ) * Y5 ( U 
1-eO-YJtJ-I Y-7-C l )-- 


c 

c 


USE 5 P-0-IN-T-- 

105 DO 111 I o j .MV 
— 1 — I— V — Y-ft-i 

T » T+H 

■E'ALL f-po - 

IF (NH) 3q 5 , ftft 9 , 6 ! 1 
IF 4-N ft - ft-) 9-5 1 r 201 . -2-1 1 
N 6 « N6+1 

. — ---N-4) — •*— N 

print 9.H,N5,N6,Nft 
9 format ( 3HH • »e t 2.5 * 5 h 
If ( N6-5 ) 500,450.450 
202 ASS1-GN 99 - TO N5W 


EOR-MUL-A -0 N G-E — S-T- frR-T-E-D 




107 

64 1- 
4 49 


NS « . I ft , 5 H M 6 * . I 5 » 5H Nft «,I5) 


C 

-C 


G0aD-P-0-f-N4- --f>GWN-54H-I-Fl" AND MOV-E'ON -- 

500 Do j20 I * J tNV 

Y~74- ! ■>— --Y 6 (’I 1 







1 1 9 * Y6"f n * ys in 

— L-2-0* Y-E44-)---.~4l%4~!~) — - = 

1 2.1 * Y’8 < I ) = Y3 ( n 

—I -22* — ; Y-3-H- ) - -«- Y - 2 -H -) 

123* Y2 < 1 ) • * Y 1 < 1> • ' 

*“1-2-%* -Y-| — (— J— ) — a — Y*0 <~H : * 

125* Yom • y ( n 

1 >~ ~ « D Y~ & ( 1 1 ■ 

1 27* DY6 ( I ) «* D Y 5 1 ll 

■ -i-2-8* — - D-Y-5 (-H * -v-CHK-H- >- r . — 

.1 29* PY8 (,I ) *» D Y 3 < I > 

.1-30* — .--DY3-U4--*-4>-Y-2-(-T-> = 

1 3 l ? * 0 Y 2 ( 1) * 0 Y.i ( J ) 

.. . 4.3-2 DXl-.t-U — * --CLlO-i-U : i 

1 3 3*. 1 20 DYO (I) « DY ( I r 

,.134^ . ‘ GO . T.O .10 5 . --- - . . , . _ ; 1 

1, 3 5 • C ' ' 

-1-3 6 • C- - NEC D- -5— P-A-S-T— F-OIN -T-5- P-L US— P-B ESE-M-T-— ROI-N-T- -F-0ft--F4-F-%-K-P E - R 1 -V A-T-i-V E- 

137* ' <150 CALL D5Y 

- -1-3-8* ; C- : - t — 

139* C D5Y RETURNS El « M A X | A B S ( Y V ) * D L I M- A B S ( Y > « E L.I M > 

1*8 0 *“ C — - - -A-NE7-E-2"*--MArXM-A B-S-(¥iM-*-DtrJ-M- A-S-S^Y4-rE-trLMy-%©-,-} 

18 1* - GO TO. NSW, (98,99) 

-1-8-2* C 

18 3* C T-HIS SEQUENCE ONLY PERFORMED IN STARTING 

— HH * ~9-8 — H= — f-Ei — g -L l -H-) -2O-2-rg-O*g-rg~0 T ' ' 

1 8 S'* C 

-1-8-6-* -C- A-L-LOVr-t'TM-'SU-CCE'SS-T-VE'-S'T-A'RTI-MG— F-/rmrRE-S 

187* 95 1 ' NF » Nf+1 

— 1-8 8* ‘ -1-F- — (N F- L-l-M-1- 9-5-2-, -9 9-l~,-9-9l r ; 

189* 201 NF a NF+1 

- - -1-50* 95^r^l-r^9--l 

1=5 1* C 

— 1_S 2-*- . C- .-.-REDUCE — g-T-E-P— AND--RE-5-T-0RE- -I-N-l-T-l-A-L—C-OND I-T4-GNS — F-O-R — N8 — * — 3 

153* 952 H * H *o , l 

-- 45 4 *- --- 00-- 9 5-3- -l*-|-f-N-y — 

-" 155* - 95 3 Yfl.( I) * Y3U ) , 

. 4-5 6-* JlO — 40 -4-9 ‘ 1 ! 

15;?* C REDUCE STEP AMD RESTORE INITIAL CONDITIONS FOR » 1 

1-5?* 950 H * H * 0,1 

-16 0* — -• -'D0--]-6G- I-= : -l ,»-N-V 

16J* DY3 ( I ) * dY8 ( I > 

. — 1-6-2-* ; Y-3-(-H— <*-- Y-8-M-) 

1 63 * 160 YO ( 1) * Y8-M ) ' 

165* ' 99 IF < E 1 * - EL I M ) 212,212,211 

167* C ERROR OUT OF ROUNDS - HALVE STEP SIZE AND TRY AGAIN 

— i~6-8 2-H — H- -a- -H - * * — 0,-5 -. — — - — 

1*69* IF (H-HMInF 990,980,980 

1*7 0 * ' - -9*80 *T- «--T- --- -5 -,-0*H'— - 

1 7 i * DO 130- 1*1 iNV 

•1*7 z*-'\ ■- * * «NAX- «-Y-l-i 1) “” : * * * * - . ' - - 

173*. _ .SET .* Y 3 ( I ,) _ - ' - - 

7 <* * Y-H-H*-H3 5*,- *’Y 0 H-T+-HC ' ♦ Q M A*X*---?©-r*Y 2-<H-*-20-,-*<j E*T— 5~, *y-4 -H--J- N*tO0*7*8 *1-2-5— 

.175* Y3( I )«(-5**Y0( n+AO»*-Q«AX + 90 . * Y 2 ( I ) -20. * S ET - ♦ 3 ♦ * Y <( ( 1 > 1 * , 007 8 1 25 

17 6 * . - -' - 


1 77* DYM { I ) = 0Y2 i I 1 

-1-7-8-* 0-Y-+4-H : 

179* 1 30 Y2 ( !) »Q«AX 

is?* caul state <Y3»dY3) 

192*' T - T+H + H- - 

183* CALL STATE < Y 1 ,dY 1 ) 

- 1 8 4'* — r - -* — -T-+-H — 

IBS* G 0 TO 10S 

187* C Error in bounds - check WHETHER TO PRINT 

1 8 8 * 212 IF t T-TPR ) --3Q2 , 30-1 I 301 . ------- 

189* 301 CALL PRYNT 

19q.*„ \ IX LX.ERj.-XF-) 21-2-i-2-i-2- t -3-Qil __ __ 

1*1* 3 0 M IF (T-TF) 302 t 306 * 306 

1’3» C ' CHECK WHETHER TO DOUBLE STEP SIZE - DELAY OF 3 PASSES 

1 9*L*. C 1-5- ~M A-N D-A-TX-D-B Y— S-T-R-U C T-U R-E — 0-F — P-04-N-T-*-S A V-* NG— ME-GH-A N+S-W 

18 &* 302 IF IN&-B) 500,680,060 

.. L-9-6-* JL8.Q.. _LF_ J-E-2-* £D.i;B.) 22-2- + .2-2.2-»-5X> 0 

197* 991 PRINT l,Nf»El 

- 1 9.8 * 1 .FORM AX (_LH £L»- J-2.,.28 H.-S-T-A R-T I-tLG— F-M VU-RE-S— «- . 

199* 990 PRINT 2*H,T 

— 2-0-0 -* — 2 .E 0-R-M AT -( 2-I-HO-S-T E P -T-0 0 - 5 M-A- L L- ■--H---«-i-E4-S »-7- t -8-H--T--«- r E -l =,-«24 

20 1* 305 PRINT 76 , CY ( I) , YOI I ) . Yt ( I > *YZ(I> , Y3( I ) ,Y*M I) ,Y5l I) .Y6U) » I«1 ,NV > 

--2-02* 7-5-f JlELMJLT4-l-FLQ-i-6.E4.7-.-8-j 

203* ■ 306 PRINT 7,NUM 

'—2 D H * 7-F-O R M-A T-<-26hG5-T-A~T-E' -X QG -ft-T-h fvN 0 -- F- V-At: ! XA-T-XD- r 1-Srr6t{— T~1 +V ERl 

2D5* PRINT 8,JAC 

216* 8 .F'O.RMA-T Li 9.HCLJ.A-C .0-8.1 J.N- E.V-A-L-UAlX£Li-t-5..»6-H T-UiE4>-) 

207* STOP 

—20 0-* 6-HO— 


ND OF COMPILATION; NO DIAGNOSTICS, 





; .33 



-SU B-R-0 U-T-I-N E~- F-P-D - 

COMMON YO,DYO»DB»T|H|NV,NUM l Y,'n ,Y2iY3»YM ,.E l' ,E2 , D Y » DY 1 .DY2.DY3, 
» .-R-T-^w 


COMMON YO ,DYO , DB fT ( H , NV ,NUM 
1 — - -O-Y -vR-T-pNO-S-X- i- U - A C Q-i~G i N- 1 

CQMMON/BNOS/YLIMiXNUIM.ELIM,! 
-D-0 UB-L-E- -PR E C 4-5-1 G N— ' f 0-4 V-O-K-O YCH- 


■ uwm, - L ■ n l sm, in IC.I 

YR(|0),OY(lO)iOYino)»DY2<lO)»DY3UO) iDYrTTo ) !T<To ) *G <1 0 * 1 0 > , 

- -E-X T-R-a-I -J-Q-) - 1 - 5-A-V-€-(-l-0-t-l-0-4-i-0-RH' l-0~)-rB£-S-H 0 I ’ i -tH-fO trD B~-dl-O-T~,C-UP-O-Y-(~l:-0"I--i 

B (to ) 


S E T -, 5 U M -r-S-t-V 

EQUIVALENCE (NV , N V 2 1 i ( Y , B ) ,(BL!M|YL1M) i(OLDB»EXtRA) 

0 A T-A / 6 N-HL ! M-/-J-7-/ — 


NB-- COUNTS E- V A LU-A-T 1-0 M 5— 0 £ • B- 

r N 10 FLAGS TWO SUCCESSIVE NOEX STEP FAILURES 

lik . X.Q-UNlS _..I±l£-5£ ...N.DO __5_T_E.E.S 

N5 “ J 


I 8*- - N 1 0 «-l - • - — — — 

1 9 * ASSIGN 05 TO NHEX 

21* ,C FIND SQRT OF PERFORMANCE INDEX, TO AVOID OVERFLOW 

- 2-2* QA-L-t— S-r-A-T-E— I-YVU-Y-) — 

23* CALL QX(DQiQMAX) 

-2H*- - -P M-»0 » • 

25* DO 20 -1*1 , N V 2 

-26-*- - ---2 0 P-I-l- » P-L l-t-(-OU <-I-I-/ QW-A-X+-*-* 2 

27* P I 1 * SORT ( P I 1 ) *0MAX - 

-2 ft* S-0— 0 0 — 1-6-3 — 

29* CURDY (IJ * OY(II 

-30-* j-6 -3—OL-DQ-H-t- B-H-) 

3 1 * 50 N6 » 1 

..0 2* — - - -,-I-O-U-T - -» — N-Y 

33* DO 25 1 * ! ,NV2 

— 3 *4.* S-U.H.3.&-,.p.Q 

35* DO 35 J * 1 , N V 2 

36* „ J35 SUM*. SUM *5A 1 ,.JJ. SJ5JUJ4 

37* IF<DA8SISUM)i*XNLJM*ISC + ABS<0LDB( 1 ) 1 J 1199, J9?,2 5 

3 ft-* ... 199 1 OU-T-* -I-eU-T.—l 

39 * 25 DB (I) -SUM 

. ‘in* .c : 

^ 1 * C IF ALL IN BOUNDS, PRESENT R«S ARE THE SOLUTION 

*2* 1.F.U OUT + -2-1-0-1-210-, -4-7- - 

*t 3 * 4 7 BN * 0 « 

HR* - — --08-W- 3--Q-* - -- 

HS* 00681*1, Nv 

R 6-* - --- - B-M-a- fi H-* — EM-l~)-»-84-I-j - — 

H 7 * 68 D8N » DBN + OB(I)*DB(I) 

- He* -- ■ - -~LF- (-DBN-BN-4- -4 3-, 4-3-,-AR 

H 9 * 43 T I a | , 

• -So* - - -- - ---60— -TO 60 ■ - — — 

5 1 9 6*4 T I * SQRT ( Bn/DBN > 

-5-?* —6R- T - 0 -6 0 1— 

. S3* 210 °0 205 I*i» NV 

• 5<4* --205 0-Y-( I ( -m CuR0-Y(-4> --- - - - ~ 

55* NH = 0 

56* - RET-UR-N — - . - — : 

57* C . 

- 5ft-* --C -- — J-F— N O-R-M-O F- - DB--T-OO--L A-R.GE-,- R-E-CA-L CU-L-A-T-E--8 — A330 — 6 F-Y— N &iy — p — I. 

5 9* 90 T I • , 1 *T ! 

40-*- . *0-0-0— 54— 1-» 1 r ^V-2--- — — - — - ~ 




-5 5 BU)>OLDB(I) - TI*DB.m 

■— — n^»4sa-j- 

• N 6 ■ N 6 + 1 


IF NDE X FAILURES, CALCULATE N 
- - ~ -Lp-(;N 6~~ NCrg-X-) — 7 0-rTOrTS 

75 IF(NIO) 230,230.220 



7,7 FORMAT ( 26H0N0 C ON VE RGENc E , N 1 0 *0 » 

--3 os- p Ri-MT— r STi^'m-rro-t , m r rt -i-rrrfr tr>~ 

7 8 FORMAT < J.X , 6 £16,8) 

■ — ' — * -LF- — -(-NH -NHL- 1-M-L---57-B - ,-S7 6-,-5 -7-6- — — 

575 NH * NH+I 

-H E4-W8H- - 

57 6 NH »' -1 

.return -- 

220 ni.qbo 


EW.G ' -• halve if THIS FA 


8H H * 


* £ I 5 » 7 
■-nTTV-f 


, 8 H T E 1 5 • 7 ) 


MUST RESTORE 

DO — 1-6! — ,-N-V 

16 1 B { I ) * OLPB { I j 


OLD -B* S 


F I ND G 

DO "l' 3 I » T ', NV 2 

X » B|,II 

EH-F — »— D'A-9'S“(~X~1 — “ — BtrMI 

IF (DIF) {6,16,15 

■ 1-6-H e-f-f-) — X~ * — {-*"00 «8 

SET » J * OD 8 

. GO--T-G— {.-8 

15 s < I > e X* I .00 1 00 

-- . S F T - ■ — }-OQ.g_,-4X 

.18 CALL STATE (Y,DY.) 

C AL-L— ® x-f-Q vG M A X-I— - 

DO 18 J « l , N V 2, 

- - —G-W- r I p* { Q-(--U-) - o q-(-U-H-*-S E T 
18 SA V E ( J , I) » G | J , I ) 

—l-3-B4-l-4-*-X _ 

CA'LL - MAT'! NV ( MV2 ,G ) 

. ..JA-C*,JA-C-+-1 


• — J-TE-R-A-T E— 0 M—T H £ — R-E-5-l-DjU A L 

DO 302 M» l »L I M 

— 1+4 5 — 

' DO 81 2 J» 1 ,NV 

-- — • SUM— »— CJ-rDO — 

DO 8 1 3 K>'I »NV 

-8-1-3 — 5-U-M — s— 5-U-M ♦--■ S A V-E-(- !■ , K- )-• ‘ 

8 l 2 R E 5 ( j , J ) a -5UM 

_,_8.|.c r .„RX-54-l-,-I-)— ■ ^-E-5-L-l-y-l-) 

- DO 820- I« j ,NV. 

- ■ 00-8-1-7- U*-] rM V : ' “ “ - 

SUM 'a 0 • DO 

D0-8T8'-K-wi V.N’V ~* ' ^ 

8 1 8 SUM » SUM + G(I,'K)*R£5(K»J) . 

— 8 -i r 7 — O &-4-J-) — -* — SyM — - 

po 81 9 *i» { »NV 

— 8 -}-9 — 5»-f-I“j-d~) — *r — QB-(-J-) — ■ — 



35. .. 


1 19* 820 CONTINUE 

-1-2 0-* 3 Q-2—C-Q NT-4 -M U£ 

121 * GO TO 50 

123* C TO HERE IF LESS THAN NOE* FAILURES 

■1 24-* 7-0-CAL-b- -S-W-E- - t-Y~,-0-Y J — — 

125* CALL QX(Q,QMAX) 

■4-2- 6 r* — P-4 -2 * g -« — — 

127* D0A0T*l,Nv2 

128*- - ‘ 8 C 1 P 12 *P I ?+( 0 ( 1 ) /©M-A-X > **'2 - 

1 29* PI 2 * SQRT < P 1 2 ) ♦ QM A X 

-130*- -C- - . 

l 3 1 * C CHECK NEW P 1 

A3.Z* LEL1E-LUEJL.U- aZjJLD».9 0 

133* - 67 GO TO NHEXt( 8 S, 86 ) 

138*-—- - 86-- IF l-N-5-200) -24Q-,-25o-r2-5Q 

135* 250^PRINT 260, H,T 

137* GO TO 305 

4- 3.8.* 8 -S—LF — LNSj»,ii£434 — 2.8 Q., . 2-8 Q-»-28-g 

139* 245 ASSIGN 86 TO NHEX 

. 48 G* GO-4-0 -7-5 

n t * c 

-1-4 2*.. C. I JTJE.R A_T £_ T-0-^fi.E-T — N ELW G-AN-0. TRY AGAIN 

143* 240 PI 1 =P 1 2 

— 1 - 8 - 8 -* N-4-O-a-l 

145* DO }10 1*1 »NV2 

- 1-8 *-*• II Q -0 ~ 0 Q { -I -)-> 

147* DJVaO»DO 

• 4 -8-8-*- - — -—GO -i-5-0 — I «-l-rN-V-2 

149* SUM»0*DD 

. -i-So-# GO— j-5rl~U*i-rN-V-2 

l 5 1 * 151 SUN = SUH + Oa< J>*Gf I , J) 

-4 52* __iUJ/-aX)4-V-*DB-n-L-*-SUM- — 

•153* 150 EXTRA!!)* SUM + T! * DB ( ! 1 

• 1 58* — — • — O-J-V — *» — «?D'!'V — — 

155* 00 155 1*1 , N V 2 

-1 -5 6-* — - — — 

157* DO t 58 J*l ,NV2 

1-58-*- 1.5 8 -SAMJE-U ,JUflflUX*SJjM. : 

159 * 155 Save u , n *1 .+save 1 1 , n 

1-6 0* D-0---4-56— K *4-»-N-V-2- 

161* 00 157 1*1, NV2 

• 4*6 2* — -S-U-M-* G»-f>0 * 

163* 00 1 58 J*1 , N V 2 

1-68* - - 1-5.8- S U M * S U H* S- A V-E- ( -H -J-)-*-G4-U-rK-} 

165* 157 EXTR A ( I ) »SUM 

— 166-* -- -DO— 1-59- -1 *1 »-NV2- 

167* 159 G( 1 ,K ) *EXTRA ( I ) 

•- 1- 6-8 •» 1-5 6 - -DQ ! )0 — ) — — 

169* GO TO 80 

1 70 * END - • ■ - - - 

:no of compilation*. no diagnostics. 



- 3 . 6 .. 


■ — J-H*- Cpr-yw-T 

2* SUBROUTINE PRYNT 

■ - - 3-* -COMMON- YO r D-Y-Ch» Ob TY-rH-r^r^It H r^rYM-T-Y^»'Y^rY^rf1'T^2-rOY-rfrY-^r9-Y'2-rO-Y^- 

r* l dyh.tpr.dt 

-e,-* — — _ — O&u BV£~-P-R^^-^ rONh ---Y ^ -<-- t - 6H -SY-6^- 1 - 0 ) - > 4 1 0 fr -Y-H-Q-H-Y- H-H I ) % -Y-24 - 4- &- W Y 3 ( - KH-y 

6* 1 Y1 ( 1 0 ) ,OY ( i 0 ) *DY 1 ( 1 0 ) * OY2 ( 1 0) >DY3 ( ) Cl ) »DYN ( 10 ) 

- • 7 * D I ■« E~N S-j ON -V-( -6--) -j I <-Y0-) 

8* EQUIVALENCE «z,db> 

lo* C LAGRANgIaN INTERPOLATION to make THE TIME come OUT nice 

1 2* 00 10 J-l ,6 

- T-3-* --- PROD- »"1t *» 

1 H * D * T 

- - - - OO— 1-2- -K-a-1 yb 

U* IF(K-J) 2 0 » 1 2 * 20 

1. 7-* 2-0-P R-00 — * — P 8-0OA.(.-T4?-R---D4-7-(-A-I-«-a-) r~ 

1 8 * 12 0 ■ 0-H 

2q* 1 0 A l » A l -H 

- -2 1* 0-0— JO— tfM-rN-V — • 

22* 30 Z(J).«V{t»*YlJ>+Vt2)*Y0(j)+V(3)*YHJ)+Vf‘t)*Y2tJ)+V(S)*Y3(J) 


-2-3'* 

2<4* 
2-5* — 


T -MH-6 1 *-Y-N“f-«J-I- 

PRINT RO.TPR 


2&* PRINT &0 , ( Z (I) , J ■ 1 » NV ) 

. 2-7-J* — 6-0 LO&MAJ — 

28* TPRaTPR*0T 

- - ■ -2 Y * RE T-U-R-N- — — — 

30* END 


3D 0J..3CLI1 PJ-UAXI.DN_ 


■ 37 - 


1* . CGILL 

- ■" 2* - — SUBRQUT-J-NE-6-I-t U — ~ “ — * s 

3* COMMON YO »OYO *U , T ,H , NV ,NUM 

~ — M-*~- .-DO U^ir£— Ffr^&l- s i - e - N ' ■ m H' O -)- >- 0 ¥ t H -tt rl ~r o t- , : 

5* Data A/»5DO,»2928932l88l3*152M76DO,l»707l06781 lS&5H752Do,*j66666666 

-- 6 • * | 6 A6 &AA b A-7-D 0-, O -2-9-2- 8-9 3-2 -1 -8-R-1-3 4B2-4-7-4R<}. 1- .- 7 -0- 7 -H3 A-^J4-S-A*-0-75-2-Ofl-t^3-3^-3 — 

7* 23333333333333 30 o,l» [)0 ».58S786*t37& 26 9o**95l DO, 3. *HM21 35623730950500 , 

B* 3 2-*-0 »-0 0 , « 2-1. 3 203^3559 A-M-25-7-3 f>Q r H t4-2-!-3-2 0-3M-35£ j ?- 6-<4-2-57-D&, O-r-Oo-/ ~ 

9* C 

1-0* :JC E-QiLEUi- I)RD-EA--BiJUG£--TaC4jJJ--A-^-eXh.0.p , RESULT- ]H Y-Q ft-RR.A-Y ..UPO-N.-E-X-U 

11* CALL STATE ( YO.DyO) 

12* - DO 500 JM»NV~ 

1 3* 500 U ( J ) *0 , DO 

- 15* 10B0 DO 1050 J«1 ,NV 

lb* Y Q-4^M— *— y-o -LJ +-A -Wf T H-*WW4--- — — — - 

17* 1 Q5 0 UJJ) ■ A ( Jl ,3) *DYO(J) - A(Jt,<n*lMJ) 

- 15*- -- • - - -GO-'T-O— 1-1-0 A-0-, 1-0-1 O-,--l-OA0-,95O-l-,-d-l 

19* 1-060 T»T *H*Q«5 

2-Q-*-- — 1-01-0- -tH — *— =J1- +--1 

21* CALL STATE! Y0,Dy0) 

6 0 — T O -l-O^'O — 

23* 950 CALL ST ATE t YO » 0 YO ) 

25* END 


l N-0- -0 F-COM fU-L.ATU.0ili JiO 04 -A-SJJ-QST-Lg-S. 


1 * 

2J* 

3* 

-* 4 * - 

5* 

7* 

8 -*- 

9* 

10*- 
1 1* 

- 1 - 2 -*- 

13* 

-1 9-* 
15* 

1 fr-*- 
17* 
-•Ja-*- 

1 9* 

'2o*‘ 

2 t * 
22 * 
23* 


COSY 


- -S U-B R.OU-T-I-N E- -0-5.Y- 


COMM.ON YO ,DYO,DB ,T,H,NV ,NUM, Y , Y1 » Y2 , Y3, Y*1 ,E1 ,E2 
— Oou&L-E ■ PR E C I-S J-0 N YQ-H-CH -j-D-yo^ T0-}-,B8-(-l'O'l - irY ClOl— j-Y-if-l-G7',-Y-2-l-l-Q-)-v'Y-3-l~l-0-T-| — 
! V H { JO) 

- 

DATA DHM/.0805/ 


-C - 

C 


30 

*1"0 

90 


20 


choose maximum truncation error (sort on- 

-E-l— »— -h,QE-20- 
£2 * E i 
_0 0 — 2-0 — U-x-i_ r t4y_ 


ER = ABS( ylvM+5,* (Y3< J)*2,*(Y! ( J 3 - Y 2 { Jl ) * YO < J ) >"«YM I J ) ) * OLIM 

E S=---A B5-(-¥-(-U-H -*-E l-l-H— - - — 


ET = Er-Es 

-l-f-(-E }-~e-T-)- 3q-,-}-0-j-1g 

El » Et 


IF{E?*ET1 Mo, 20,20 

E2- Er 

CONT l NUE 

RETURN - - ... 

END 


;No -OF COM-p-I-L-A-TlON-: No DJ.a 6 N OS-T- 1 C S-, - 



38 


- ■ • - - C : Q X-- * — 

2* - SUBROUT I NE Q X ( 0 | Z ) 

-3-* C-Q B -M O - N ■ . Ve -T-a-Y-O-r ^ -r^-ti-rN-V- t^l - u - M ' -i - V -, Y l -,- Y2 , Y3 , Y R , E 1 , E 2 - ,'P Y r frV i r- P - Y2 - v fr Y- 3 

R* DOUBLE PRECISION YO t 1 0 ) , D YO U 0 ) » DB M 0 ) i Y M 0 ) i Y U 1 0 ) » Y2 ( 1 0 ) » Y3 < 1 0 ) » 

.km - 1 ¥H-L J 43 ^ r 0^t^0 ^M-I-MD-Lt-D-Y-2 <-I-04-i-0 J ( , -3-H-G-> 

6 * DOUBLE- PRECISI ON Q,ZiX 

- -7*~ D-I-M-e-NS-J-Otl-^-H-I 

8 * DOUBLE- PRECISION BH,BO,B2,B3 t Bl .AO.A1 ,A2,A3 

< 5 * P A4- A. --A OT -Kr-S-B-%2-^9 9 -9- 9 -R^-9-9^-6^K>-0-y ? 

ID* DATA Aj /-l .Ot 696529999996*71 Do/ 

~‘M* DAT-A A 2 / # 5 - 2 9 R R-R 0 R 9 9 9 9 97 0 9 8 8 D 0 / 

12* DATA A3/- ,-967787R9999?9 1OR630- ! / 

■13 * O-A 7-A- -B MV ."R^-3-B-6-6“R-2-2-H ?-*» 7-J-3 R-SiJOy-- -• 

1*(* DATA Bo/ • 235026309RR67B72R3D0/ 

.... 15 , DA4A-6i-/ T ^'6- W 7' fr < t »< P 2 ~ 6f3 2 R fr^-fr/ 

1 6 * DATA Bz/. 567200, 5 28OJ0787R05D-2/ 

--- 1-7'*-“ DAT-A - B-3-/-.-2-0-1 3 3-67-0 0 875 9-R-7-6 -2-7 2D-3V 

18* C 

- -19 * C - • T-AKE D-I FF-E-R E-NC-E-0 F -Y--R-R E D-I-C-T- E-P-4-0 R4-S-J-N-A L U-T--f.ROH--EXXRAP0L-A-T-I-0.N4 — 

2q* C AND-Y CALCULATED ( 8 A 5 EP ON INTEGRATION FORMULA) 

2-j* Z-— *— G-.-B -0 

22* . DO 30 U*1 ,Nv 

2 - 3 *- — Q { J ..y-(-j-)-*4-y-0LJ-)-*-A 0 + Y- ! (J-)-*A i.+.y. 2 .( J-)-*-A- 2 + Y-34 U-!-*-A 3 - 

2r* 1 H*(DYI J)*BM + 0Y0( J)*BO+DYU J)*B1 +DY2I J)*B2 + DY3( J)*B3) ) 

-- -25* — X--=-D-ArB5-(-Q-0U-)-) - - - 

26 * IF < X - Z ) 30 ( 30 1 R 0 

2 7-* -R-0- 2—-* — -X — 

28* 30 CONTINUE 

2 9 .* 4j_EJ4ARjs) 


1* CMATINV 

“ - 2* - • ■ - SUB-pOUri NE~ M-ATi-WV-(N-iZl * ' -- 

3* DOUBLE PRECISION Z.X.CC 

i,-* -D-+-M E-NS4-0N-“hp-H--0-nZ-<'t-8-rt-r-fi-X-<--ltt'l“2-) 

5* DO JO Is) t N 

6*' ~ tO”I P<*t-)'*0" ' ‘ “ 

7* DO 76 1*1 f N 

... fl .» — X — *— 0-r00 

9* DO 69 J= I ,N 

1 1* 6 ! DO 68 K*1 ,N 

• 1-2-* - IF t-IP-lK- 1--T-) -6R-,-6 8-,-3R9- 

13* 64 CC»DA8S ( Z { J ♦ K ) > 

1-4* I-F I X-— -0C-) — 6-7-»-6 8- r 6-8 

15* 67 1 2*U 

1-6* — — • — 

17* X * CC 

.. .18*. . - 68 X QHX.LNUE. 

19* 69 CONTINUE 

-2 0 *- - i-PR Pi -H-4-+-H - 

2 1 * I F I 12-1 1 ) 7 0 » 7 2 » 7 0 

— 2 2* ~ — r 7-0~D-0“"7i — 

23* X«Z ( I 2 , L 1 

24* ■" - z-( iz *Li*?rTr,L ) •--- — 

25* 71 Z « I J » L ) * X 

26* - -7 2' "I X ( I » 1 > *T2 * ■ * — • 

27* I X ( I , 2 ) * I J 

- - — 2-3* — X « — -Z — (r-t~ ) — f — l } — — — - 

29* - 1 ( I ! * I 1 ) *1 .0 

~ 30* DO — 7 - 3 - )-j-N - • - - — 



3 1 « 

- — 3-2 * - 

33* 
— -3 4 
35* 

-- 3 6 * 
37* 

-3*}-*- 

39* 

Mq*., 

H]* 

-‘♦•'2'* 
<4 3 • 
.. „ iU|.* . 
*45* 

. - 4 6 . 

. H 7 * 

-^e* 

H9* 

So* 

5 l * 

- 52* 


73 Z|I|,L)«Zt!liL)/X 

po— 7-6- -J " 1 r N ■ — 

IF < J-I 1 ) 74 ,76 ,7*4 

- - • 74-Xr«-Z~<-J- r H--) 

zu.it > * o * 

DO .75 — lr*-l-»N 

75 Z( J,L)=Z< J>L>-Z< It »L)*X 

76—C^-f-i-T-J-N U £ 

L * N+ j 

DO- 79- i -4 t N - 

L * L- 1 

— - ■ I F. < !-X (-t. »-!-> — I X -7-7 j7 9- r 7-7 

77 I 2 - I X t L , l ) 

- — l~l~»4-)U-L~t-2^ 

DO 78 K* 1 » N 

X*Z.<K,124 - - - - 

Z(K,12)«Z{K t !]) 

Z.t.K.^44 )*-X- 

78 CONTINUE 

— 7-9- CotlT-I-N-U-E 

399 RETURN 


•ND OF COMPILATION 



NO diagnostic 5 . 



C ST A T E 


COMMON VO » DyO , D0 » T iH »NV ,NUM 
— Q-O-U BT.-E — P-H-E-C-I-S-I-ON — T6-f - t~Q-)-r^T~9 ( 1 0 T~vO B ( 1- { 
DOUBLE PRECISION X< 1 > ,D( 1 ) 


pO{ 2)-*-XU-)*X (2T 
C „D(.l )»-lOOO.*( X ( 


0 ( 2 -)) 


--num-_»--n-u-h-+-: 

RETURN 

end 


OF COMPILATION; 


DIAGNOSTICS* 




0 . 1 0 n ;iO(i o 01 (! . 100D0006 01 


kl 


LIMITS' 

' 0 . 100fl0n0i-05 

H = O.IOPQOE-Ol 

h = n.iononfc-oi 
h = o , ionont -02 
h = o.ionnoh-o 2 

H = 0. In 00 0b-03 
H = 0 . 100006-03 

m = o.ionortfc-o3 

H = 0.10000E-Q3 

M = 0,100006-03 
» = o.ionoob-o3 
h = o.20nont:-o3 
H = 0 , 200006-03 
H = 0 . 2<l0.0 0E-0 3 
H = 0.200006-03 
H = 0.200006-03 
H= 0.200006-03 
H‘ = 0. 200006-03 
H = 0.200006-03 
H 2 0.200006-03 
H = 0.200006-03 
H = 0.400006-03 
H = 0.200006-03 
H = 0. 200006-03 
H = 0.200006-03 
H = 0, 20O00E-03 
N = 0.200006-03 
H = 0.200006-03 
H = 0.400006-03 
H 2 0.400006-03 
H = 0,400006-03 
H = 0.400006-03 
H = 0.400006-03 
H = 0,800006-03 
H = 0 , 800006-03 
H = 0.800006-03 
H = 0,800006-03 
H 2 0,800006-03 
H = 0.160006-02 
H = 0. 160006-02 

w = 0.160006-02 
H = 0. 160006-02 
H = 0 . 320006-02 
H = 0.320006-02 
H = 0.32n0 06-02 
H - 0.320006-02 
H 2 0 . 64nOOE-0? 

H 2 0, 640Q 06-0? 

H = 0 . 640006-02 
H 2 0.640006-02 
H n 0,640006-02 
H = 0. 640006-02 
H = 0.12P006-01 
H = 0.12fl0»6-01 
H = 0.1?fi0 06-01 
H = 0.1280 0 6- 01 
H = 0.266006-01 
H = 0.254006-01 
H. = 0 . 26 61106 - 0 1 
h = 0.2^61)06 -01 


0 . 1000OQ 0E-03 
N5 a 26 N 6 » 

N5 2 2 N 6 » 

N5 * 2 N 6 » 

N5 2 2 N 6 2 

N5 = 2 N6 2 

N5 2 2 N 6 * 

N5 2 2 M 6 2 

N5 2 2 M6 2 

M5 2 ? N6 » 

N5 2 2 N 6 = 

N5 2 3 N 6 2 

N5 2 2 N 6 = 

N5 = 2 N 6 = 

N5 2 2 N6 = 

N5 3 2 N 6 = 

N5 2 ‘ 2 N 6 2 

N5 =. 2 N 6 = 

N5 = 2 N 6 2 

N5‘ 2 2 N 6 = 

N5 2 2 M 6 2 

N5 2 3 N 6 = 

N5 = 3 N 6 = 

N5 = 2 N 6 2 

N5 = 2 \i 6 2 

N5 = 2 N 6 2 

N5 « 2 N 6 = 

N5 2 2 N 6 = 

N5. *• 3 N 6 = 

N5 3 2 N 6 2 

N5 3 2 N 6 2 

N5 = 2 N 6 2 

N5 2 2 N 6 = 

N5 3 3 (j6 = 

N5 2 ? N 6 2 

N5 2 1 M 6 s 

N5 2 2 N 6 a 

N5 2 2 N 6 = 

N5 3 3 N 6 3 

N5 2 2 N 6 = 

N5 = 1 N 6 - 

N5 2 1 N 6 = 

N5 2 3 N 6 2 

N5 = 1 N 6 = 

N5 2 1 N 6 ® 

N5 = 2 N 6 2 

N5 2 3 N 6 * 

N5 = 1 N 6 = 

M5 = 2 N6 = 

N5 2 2 N 6 ' = 

N5 2 2 M 6 = 

N5 = 2 N 6 = 

M5 = 3 N 6 = 

N5 = 2 -M 6 = 

M5 2 2 M6 = 

^5 = 1 N 6 : 

N5 = 2 N 6 = 

N5 = 1 M 6 = 

^5 2 2 M6 = 

N!<3 = ? M6 = 


0 . 1 0000006-04 

4 N'3 4 

5 N4 = 5 

4 N4 « 4 

5 N4 « 5 

4 N4 ■ 4 

5 M4 « 5 

6 N4 • 6 

7 N4 ■ 7 

8 N4 • 8 

9 N4 * -9 

5 «4 « 10 

6 N4 ® 11 

7 N4 « 12 

8 M4 • 13 

9 N4 2 14 

10 N4 *' 15 

11 N4 «? 36 

12 N4 2 17 

13 N4 « 18 

14 N4 2 19 

5 N4 3 20 

5 N4 ■ 21 

6 N4 1 22 

7 N4 * 23 

8 N4 2 24 

9‘ N4 a 25 

10 N4 a 26 

5 N4 * 27 

6 N4 • 28 

7 N4 ■ 29 

8 N4 2 30 

9 N4 a 3i 

5 N4 « 32 

6 N4 » 33 

7 N4 2 34 

8 N4 « 35 

9 N4 2 36 

5 N4 * 37 

6 N4 " 38 

7 N4 2 39 

8 N4 2 40 , 

5 N4 " 41 

6 M « 42 

7 N4 ». 43 

8 N4 2 44 

5 N 4' 2 4 5 

6 N4 2 46 

7 N4 b 47 

8 N4 2 48 

9 N4 ■ 49 

10 M 2 50 

5 N4 a 51 

6 N4 s 52 

7 M 3 53 

a M * 54 

5 N4 « 65 

6 N4 = 66 

/ W4 * „ 5 7 


U, 100000 HE 01 



3 



x 

n . 256«oe-oi ^5 = 

2 N6 = 

9 

N4 

5 

59 

T = 0. 2500000b Cl 





H — 

n.25*,noE-oi M5 = 

2 NI6 = 

10 

N4 

= 

60 

0 .143365 

0,286662 




1-i 3 

o.236ont-or ^5 = 

1 N6 = 

11 

N4 

s 

61 

H = 0.40960b 00 ^5 = 

4 M6 = 5 


3 

6tJ 

ri — 

0. 51500 E* 01 ^5 X 

? M6 = 

5 

IM4 

X 

62 






H = 

0 . 51 ?ont-oi ^5 = 

2 N6 = 

6 

N4 

.5 

63 

T x 0.300UD00E 01 





H x 

o. 5i?nofc-ni N5 = 

2 N6 = 

7 

N4 

X 

64 

0 .111655 

0,223256 




H x 

0.51501)6*01 N<5 = 

2 N6 = 

8 

N4 

3 

65 

H x 0.40960b 0-0 ^5 = 

2 N6 = 6 


= 

61 








W = 0.40960E 00 N5 = 

? W6 = 7 

M 

X 

85 

T 

= 0. 5000000^ 00 












0. 389622 

0 . 779044 





T = 0.3500000b 01 





>A s 

0.515006-01 ^5 = 

? N6 a 

9 

N4 

0 

66 ■ 

0 . 86 9545E-0 1 

0.173853 




H = 

0.512006*01 N5 = 

1 N6 = 

10 

N4 

X 

67 

H = 0.409606 00 ^5 = 

2 N6 = 8 

N4 

e 

6 Vi 

W = 

0.1 0740 E 00 N5 = 

3 n6 = 

5 

N4 

6 

68 






H - 

' 0 . 105406 00 N5 = 

1 M6 = 

6 

N4 

9 

69 

T = 0.4000000b 01 





M X 

0 . 1 0 ? 4 0 E 00 N5 = 

1 N6 = 

7 

N4 

9 

70 

0.677053E-01 

0.135392 




M = 

0 . 10 24 0 6 00 N5 = 

2 M6 = 

8 

N4 

X 

71 

H = 0,409606 00 N5 = 

2 N6 = 9 

N4 

X 

84 

T 

= O.IOOOOOOE 01 






T = 0.4500000b 01 






0 .303454 

0.606760 





0.5274306*01 

0.105446 




H • = 

O.2O40OE CTO ^5 = 

3 N6 = 

5 

N 4 

-C 

72 

H x 0.40960b 00 N5 = 

2 N6 x jo 

N4 

c 

b5 

H X 

0.20480E 00 M5 = 

2 N6 = 

6 

N4 

X 

73 






H = 

0.20480E 00 N5 = 

2 N6 = 

7 

N4 

X 

74 

T * 0. 5000000b 01 












0. 4106776*01 0 

.821201E-01 




T ' 

= 0 . 1 500000b 01 






H = 0.409606 00 N5 = 

2 N6 = 11 

N4 

X 

86 


0 . 236346 

0 .472573 










H = 

0.2 0 4.8 0E 0 0 N 5 = 

2 N6 s 

8 

N4 

X 

75 

. T * 0. 5500000b 01 





h = 

0,204806 00 N5 = 

3 N6 s 

9 

N4 

3 

76 

0,3198466-01 0 

, 639539F-01 











H = 0.40960E 00 N5 x 

2 N6 c 12 

N4 

3 

87 

T : 

= 0.2000000b 01 






H = 0.40960E 00 N5 ' = 

2 N6 = 13 

N4 

X 

88 


0 . 194070 

0,368061 










H s 

0.20480E 00 n 5 = 

2 N6 = 

10 

N4 

* , 

77 

T s 0,6000000b 01 





H = 

0.204806 00 N5 a 

2 N6 = 

11 

N4 

X 

78 

0.2490936-01 0 

.498066E-01 




H = 

0.20480E 00 N5 = 

2 N6 = 

12 

N4 

X 

79 







-fT 

r o 


U3 


H S 0,819206 DO N5 3 5 N6 = 5 M » 89 

T = 0.650000Q6 01 

0.1940066-01 G.387900E-01 

T = 0.70000006 01 

0.1510006-01 0.301892E-01 

H s 0,819206 00 N5 5 4 N6 = 6 N4 * 90 

T = 0.75II0000E 01 

• • 0.11 73096-01 0.2345686-Di 

H s 0.819206 00 N? * 2 N6 = 7 N4 ■ 91 

T * 0.«0O0000 £ 01 

0,9116136-02 0.182281E-01 

T = 0.85000006 01 

0.7085/66-02 0.14J683E-01 

H s D.81920E 00 N5 = 2 N6 = 8 N4 » 92 

T = 0.9000000^ 01 

0.550849E-02 0.110139E-01 

T = 0.9500000E 01 

0.4286156-02 0,8569836-02 

M s 0,819206 00 N5 = 2 N6 = 9 N4 ■ 93 

T = 0. 1000000E 02 ' 

0.3335206-02 0.666902E-02 

state equations evaluated ?/a times 

JACOBIAN EVALUATED 4 TIMES 



